raster_2.3-0_(2014-9-5_by_Robert-J.-Hijmans-[cre,-aut],)

# R classes for raster (grid) type spatial data
# Robert J. Hijmans, r.hijmans@gmail.com
# November 2008
# Version 1.0
# Licence GPL v3
setClass('Extent',
        representation (
                xmin = 'numeric',
                xmax = 'numeric',
                ymin = 'numeric',
                ymax = 'numeric'
        ),      
        prototype (     
                xmin = 0,
                xmax = 1,
                ymin = 0,
                ymax = 1
        ),
        validity = function(object)     {
                c1 <- (object@xmin <= object@xmax)
                if (!c1) { stop('invalid extent: xmin >= xmax') }
                c2 <- (object@ymin <= object@ymax)
                if (!c2) { stop('invalid extent: ymin >= ymax') }
                v <- c(object@xmin, object@xmax, object@ymin, object@ymax)
                c3 <- all(!is.infinite(v))
                if (!c3) { stop('invalid extent: infinite value') }             
                return(c1 & c2 & c3)
        }
)
setClass('.Rotation',
        representation (
                geotrans = 'numeric',
                transfun = 'function'
        )
)
setClass ('BasicRaster',
        representation (
                title = 'character',
                extent = 'Extent',
                rotated = 'logical',
                rotation = '.Rotation',
                ncols ='integer',
                nrows ='integer',
                crs = 'CRS',
                history = 'list',
                #meta = 'list',
                z = 'list'
        ),
        prototype (     
                rotated = FALSE,
                ncols= as.integer(1),
                nrows= as.integer(1),
                crs = CRS(),
                history = list(),
                #meta = list(),
                z = list()
        ),
        validity = function(object) {
                validObject(extent(object))
                c1 <- (object@ncols > 0)
                if (!c1) { stop('ncols < 1') }
                c2 <- (object@nrows > 0)
                if (!c2) { stop('nrows < 1') }          
                return(c1 & c2)
        }
)
setClass ('Raster', contains = c('BasicRaster', 'VIRTUAL') )

setClass('.RasterFile', 
        representation (
                name ='character',
                datanotation='character',
                byteorder ='character',
                nodatavalue ='numeric', # on disk, in ram it is NA
                NAchanged ='logical',
                nbands ='integer',
                bandorder ='character',
                offset='integer',
                toptobottom='logical',
                blockrows='integer',
                blockcols='integer',
                driver ='character',
                open = 'logical'
                ),
        prototype (     
            name = '',
                datanotation='FLT4S',
                byteorder = .Platform$endian,
                nodatavalue = -Inf,
                NAchanged = FALSE,
                nbands = as.integer(1),
                bandorder = 'BIL',
                offset = as.integer(0),
                toptobottom = TRUE,
                blockrows = as.integer(0),
                blockcols= as.integer(0),
                driver = '', 
                open = FALSE
        ),
        validity = function(object) {
                c1 <- datanotation %in% c('LOG1S', 'INT1S', 'INT2S', 'INT4S', 'INT1U', 'INT2U', 'FLT4S', 'FLT8S')
                return(c1)
        }
)
setClass('.SingleLayerData', 
        representation (
                values='vector', 
                offset='numeric',
                gain='numeric',

                inmemory='logical',
                fromdisk='logical',

                isfactor = 'logical',
                attributes = 'list',

                haveminmax = 'logical',
                min = 'vector',
                max = 'vector',
                band = 'integer',
                unit = 'character',
                names = 'vector'
                ),
        prototype (     
                values=vector(),
                offset=0,
                gain=1,

                inmemory=FALSE,
                fromdisk=FALSE,
                isfactor = FALSE,
                attributes = list(),

                haveminmax = FALSE,
                min = c(Inf),
                max = c(-Inf),
                band = as.integer(1),
                unit = '',
                names=c()

        ),      
        validity = function(object) {
        }
)
setClass ('.RasterLegend',
        representation (
                type = 'character',
                values = 'vector',
                color = 'vector',
                names = 'vector',
                colortable = 'vector'
                ),
        prototype (
                )
        )


setClass ('RasterLayer',
        contains = 'Raster',
        representation (
                file = '.RasterFile',
                data = '.SingleLayerData',
                legend = '.RasterLegend'
        )
)
setClass('.MultipleRasterData', 
        representation (
                values='matrix', 
                offset='numeric',
                gain='numeric',
                inmemory='logical',
                fromdisk='logical',
                nlayers='integer',
                dropped = 'vector',
                isfactor = 'logical',
                attributes = 'list',
                haveminmax = 'logical',
                min = 'vector',
                max = 'vector',
                unit = 'vector',
                names= 'vector'

                ),
        prototype (     
                values=matrix(NA,0,0),
                offset=0,
                gain=1,
                #indices =vector(mode='numeric'),
                inmemory=FALSE,
                fromdisk=FALSE,
                nlayers=as.integer(0),
                dropped=NULL,
                isfactor = FALSE,
                attributes = list(),
                haveminmax = FALSE,
                min = c(Inf),
                max = c(-Inf),
                unit = c(''),
                names = c('')
        ),      
        validity = function(object) {
        }
)
setClass ('RasterBrick',
        contains = 'Raster',
        representation (
                file = '.RasterFile',
                data = '.MultipleRasterData',
                legend = '.RasterLegend'
        )
)


setClass ('RasterStack',
        contains = 'Raster',
        representation (
            filename ='character',
                layers ='list'
                ),
        prototype (
                filename='',
                layers = list()
                ),
        validity = function(object) {
                if (length(object@layers) > 1) {
                        cond <- compareRaster(object@layers, extent=TRUE, rowcol=TRUE, tolerance=0.05, stopiffalse=FALSE, showwarning=FALSE) 
                } else {
                        cond <- TRUE
                }
                return(cond)
        }
)
setClassUnion(RasterStackBrick, c(RasterStack, RasterBrick))
setClass ('RasterLayerSparse',
        contains = 'RasterLayer',
        representation (
                index = 'vector'
        ),
        prototype (
                index = vector(mode='numeric')
        )
)       
setClass ('.RasterBrickSparse',
        contains = 'RasterBrick',
        representation (
                index = 'vector'
        ),
        prototype (
                index = vector(mode='numeric')
        )
)       
setClass ('.RasterQuad',
        contains = 'Raster',
        representation (
            filename ='character',
                bricks ='list'
                ),
        prototype (
                filename='',
                bricks = list()
                ),
        validity = function(object) {
                if (length(object@bricks) > 1) {
                        test <- compareRaster(object@bricks, extent=TRUE, rowcol=TRUE, tolerance=0.05, stopiffalse=FALSE, showwarning=FALSE) 
                } else {
                        test <- TRUE
                }
                return(test)
        }
)
#setClassUnion(RasterStackBrickList, c(RasterStack, RasterBrick, RasterList))
setClass ('.RasterList',
        contains = 'list',
        representation (),
        prototype (),
        validity = function(object) {
                s <- sapply(object, function(x) inherits(x, 'Raster'))
                return( sum(s) == length(s))
        }
)

# Author: Robert J. Hijmans
# Date : June 2008
# Version 0.9
# Licence GPL v3
.addFiles <- function(x, rasterfiles, bands=rep(1, length(rasterfiles))) {
        if (length(bands) == 1) {
                bands=rep(bands, length(rasterfiles))
        } 
        rasters <- list()
        for (i in 1:length(rasterfiles)) { 
                if (bands[[i]] < 1) {
                        r <- raster(rasterfiles[[i]], band=1)
                        rasters <- c(rasters, r)
                        if (nbands(r) > 1) {
                                for (j in 2:nbands(r)) {
                                        r <- raster(rasterfiles[[i]], band=j)
                                        rasters <- c(rasters, r)
                                }
                        }
                } else {
                        rasters <- c(rasters, raster(rasterfiles[[i]], FALSE, band=bands[[i]]))
                }
        }       
        x <- addLayer(x, rasters) 
        return(x)
}

# Author: Robert J. Hijmans, r.hijmans@gmail.com
# Date : September 2009
# Version 0.9
# Licence GPL v3

if (!isGeneric(addLayer)) {
        setGeneric(addLayer, function(x, ...)
                standardGeneric(addLayer))
}       
setMethod('addLayer', signature(x='Raster'), 
function(x, ...) {
        rasters <- .makeRasterList(...)
        if (! inherits(x, 'RasterStack')) {
                x <- stack(x)
        }
        if (length(rasters)==0) { 
                return(x) 
        }
        if (nlayers(x) > 0) {
                compareRaster(c(x, rasters))
        } else if (length(rasters) > 1) {
                compareRaster(rasters)
        }

        vals <- sapply(rasters, hasValues) 
        if (sum(vals) == 0 &  nlayers(x) == 0) { 
                vals[1] <- TRUE 
        }
        if (sum(vals) != length(vals)) { 
                warning('Cannot add a RasterLayer with no associated data in memory or on disk to a RasterStack')
        }
        rasters <- rasters[vals]

        if (nlayers(x) == 0) {
                r <- rasters[[1]]
                x@nrows <- r@nrows
                x@ncols <- r@ncols
                x@extent <- r@extent
                x@crs <- r@crs
                if (rotated(r)) {
                        x@rotated = r@rotated
                        x@rotation = r@rotation
                }
                nl <- 1
                x@layers[nl] <- r 
                rasters <- rasters[-1]
                if (length(rasters)==0) { return(x) }
        }
        x@layers <- c(x@layers, rasters)
        names(x) <- sapply(x@layers, names)
        return(x)
}       
)

# Author: Jacob van Etten jacobvanetten@yahoo.com
# Date :  January 2009
# Version 0.9
# Licence GPL v3
.cs <- function(a,b) {
        aRep <- rep(a,times=length(b))
        cbind(aRep,as.integer(aRep+rep(b,each=length(a))),deparse.level=0)
}
.adjacency <- function(x, ...) {
        warning('function adjaceny is obsolete and will be removed from the raster package.\nUse function adjacent in stead')
        dots <- list(...)
        fromCells <- dots$fromCells
        toCells <- dots$toCells
        directions <- dots$directions 

        if (is.character(directions)) { directions <- tolower(directions) }
        stopifnot(directions %in% c(4,8,16) | directions=='bishop')
        x <- raster(x)
        outerMeridianConnect <- .isGlobalLonLat(x)

        if (directions==bishop) { 
                return(.adjBishop(x, fromCells, toCells, outerMeridianConnect)) 
        }

        nCols <- ncol(x)
        nCells <- ncell(x)

        left <- seq(nCols+1,(nCells-2*nCols+1),by=nCols) 
        right <- seq(2*nCols,nCells-nCols,by=nCols)
        upper <- 2:(nCols-1)
        lower <- seq((nCells-nCols+2),(nCells-1),by=1)
        upperleft <- 1
        upperright <- nCols
        lowerleft <- nCells-nCols+1
        lowerright <- nCells
        fromCellsCore <- as.integer(setdiff(fromCells,(c(left,right,upper,lower,upperleft,upperright,lowerleft,lowerright))))
        fromCellsUpper <- as.integer(intersect(fromCells,upper))
        fromCellsLower <- as.integer(intersect(fromCells,lower))
        fromCellsLeft <- as.integer(intersect(fromCells,left))
        fromCellsRight <- as.integer(intersect(fromCells,right))
        fromCellUpperleft <- as.integer(intersect(fromCells,upperleft))
        fromCellUpperright <- as.integer(intersect(fromCells,upperright))
        fromCellLowerleft <- as.integer(intersect(fromCells,lowerleft))
        fromCellLowerright <- as.integer(intersect(fromCells,lowerright))
        rook <- c(1,-1,nCols,-nCols)
        coreFromToRook <- .cs(fromCellsCore,rook)
        upperFromToRook <- .cs(fromCellsUpper,rook[1:3])
        lowerFromToRook <- .cs(fromCellsLower,rook[c(1,2,4)])
        leftFromToRook <- .cs(fromCellsLeft,rook[c(1,3,4)])
        rightFromToRook <- .cs(fromCellsRight,rook[2:4])
        upperleftFromToRook <- .cs(fromCellUpperleft,rook[c(1,3)])
        upperrightFromToRook <- .cs(fromCellUpperright,rook[2:3])
        lowerleftFromToRook <- .cs(fromCellLowerleft,rook[c(1,4)])
        lowerrightFromToRook <- .cs(fromCellLowerright,rook[c(2,4)])
        fromto1 <- rbind(coreFromToRook,upperFromToRook,lowerFromToRook,leftFromToRook,rightFromToRook,upperleftFromToRook,upperrightFromToRook,lowerleftFromToRook,lowerrightFromToRook)

        if (outerMeridianConnect) {
                meridianFromLeft <- rbind(
                        cbind(fromCellsLeft,as.integer(fromCellsLeft+nCols-1)),
                        cbind(fromCellUpperleft,as.integer(fromCellUpperleft+nCols-1)),
                        cbind(fromCellLowerleft,as.integer(fromCellLowerleft+nCols-1))
                        )
                meridianFromRight <- rbind(
                        cbind(fromCellsRight,as.integer(fromCellsRight-nCols+1)),
                        cbind(fromCellUpperright,as.integer(fromCellUpperright-nCols+1)),
                        cbind(fromCellLowerright,as.integer(fromCellLowerright-nCols+1))
                        )
                fromto1 <- rbind(fromto1,meridianFromLeft,meridianFromRight)
        }
        fromto <- subset(fromto1,fromto1[,2] %in% toCells)
        if (directions > 4)     {
                bishop <- as.integer(c(-nCols-1, -nCols+1, nCols-1,+nCols+1))

                coreFromToBishop <- .cs(fromCellsCore,bishop)
                upperFromToBishop <- .cs(fromCellsUpper,bishop[3:4])
                lowerFromToBishop <- .cs(fromCellsLower,bishop[1:2])
                leftFromToBishop <- .cs(fromCellsLeft,bishop[c(2,4)])
                rightFromToBishop <- .cs(fromCellsRight,bishop[c(1,3)])
                upperleftFromToBishop <- .cs(fromCellUpperleft,bishop[4])
                upperrightFromToBishop <- .cs(fromCellUpperright,bishop[3])
                lowerleftFromToBishop <- .cs(fromCellLowerleft,bishop[2])
                lowerrightFromToBishop <- .cs(fromCellLowerright,bishop[1])
                fromto2 <- rbind(coreFromToBishop,upperFromToBishop,lowerFromToBishop,leftFromToBishop,rightFromToBishop,upperleftFromToBishop,upperrightFromToBishop,lowerleftFromToBishop,lowerrightFromToBishop)

                if (outerMeridianConnect) {
                        meridianFromLeft <- rbind(
                                .cs(fromCellsLeft,c(2*nCols-1,-1)),
                                cbind(fromCellUpperleft,as.integer(fromCellUpperleft+2*nCols-1)),
                                cbind(fromCellLowerleft,as.integer(fromCellLowerleft-1))
                                ) 
                        meridianFromRight <- rbind(
                                cbind(rep(fromCellsRight,times=2),as.integer(c(fromCellsRight-2*nCols+1,fromCellsRight+1))),
                                cbind(fromCellUpperright,as.integer(fromCellUpperright+1)),
                                cbind(fromCellLowerright,as.integer(fromCellLowerright-2*nCols+1))
                                )
                        fromto2 <- rbind(fromto2,meridianFromLeft,meridianFromRight)
                }
                fromto2 <- subset(fromto2,fromto2[,2] %in% toCells)
                fromto <- rbind(fromto,fromto2)
        }
        if (directions > 8) {
                leftOuter <- seq(2*nCols+1,nCells-3*nCols+1,by=nCols) 
                rightOuter <- seq(3*nCols,nCells-2*nCols,by=nCols)
                upperOuter <- seq(3,nCols-2,by=1)
                lowerOuter <- seq(nCells-nCols+3,nCells-2,by=1)
                upperleftUnder <- nCols+1
                upperrightLeft <- nCols-1
                lowerleftUp <- nCells-2*nCols+1
                lowerrightUp <- nCells-nCols            
                upperleftRight <- 2
                upperrightUnder <- 2*nCols
                lowerleftRight <- nCells-nCols+2
                lowerrightLeft <- nCells-1
                leftInner <- seq(2*nCols+2,(nCells-3*nCols+2),by=nCols) 
                rightInner <- seq(3*nCols-1,nCells-2*nCols-1,by=nCols)
                upperInner <- seq(nCols+3,2*nCols-2,by=1)
                lowerInner <- seq(nCells-2*nCols+3,nCells-nCols-2,by=1)
                upperleftInner <- nCols+2
                upperrightInner <- 2*nCols-1
                lowerleftInner <- nCells-2*nCols+2
                lowerrightInner <- nCells-nCols-1
                fromCellsCoreInner <- setdiff(fromCells,(c(leftOuter,rightOuter,upperOuter,lowerOuter,upperleft,upperright,lowerleft,lowerright, upperleftUnder, upperrightLeft, lowerleftUp, lowerrightUp, upperleftRight, upperrightUnder, lowerleftRight, lowerrightLeft, leftInner, rightInner, upperInner, lowerInner, upperleftInner, upperrightInner, lowerleftInner, lowerrightInner))) 

                fromCellsUpperInner <- as.integer(intersect(fromCells,upperInner))
                fromCellsLowerInner <- as.integer(intersect(fromCells,lowerInner))
                fromCellsLeftInner <- as.integer(intersect(fromCells,leftInner))
                fromCellsRightInner <- as.integer(intersect(fromCells,rightInner))
                fromCellUpperleftInner <- as.integer(intersect(fromCells,upperleftInner))
                fromCellUpperrightInner <- as.integer(intersect(fromCells,upperrightInner))
                fromCellLowerleftInner <- as.integer(intersect(fromCells,lowerleftInner))
                fromCellLowerrightInner <- as.integer(intersect(fromCells,lowerrightInner))     
                fromCellsLeftOuter <- as.integer(intersect(fromCells,leftOuter))
                fromCellsRightOuter <- as.integer(intersect(fromCells,rightOuter))
                fromCellsUpperOuter <- as.integer(intersect(fromCells,upperOuter))
                fromCellsLowerOuter <- as.integer(intersect(fromCells,lowerOuter))
                fromCellUpperleftUnder <- as.integer(intersect(fromCells,upperleftUnder))
                fromCellUpperrightLeft <- as.integer(intersect(fromCells,upperrightLeft))
                fromCellLowerleftUp <- as.integer(intersect(fromCells,lowerleftUp))
                fromCellLowerrightUp <- as.integer(intersect(fromCells,lowerrightUp))
                fromCellUpperleftRight <- as.integer(intersect(fromCells,upperleftRight))
                fromCellUpperrightUnder <- as.integer(intersect(fromCells,upperrightUnder))
                fromCellLowerleftRight <- as.integer(intersect(fromCells,lowerleftRight))
                fromCellLowerrightLeft <- as.integer(intersect(fromCells,lowerrightLeft))
                knight <- c(-2*nCols-1, -2*nCols+1, -nCols-2, -nCols+2, nCols-2, nCols+2, 2*nCols-1, 2*nCols+1) 

                coreInnerFromToKnight <- .cs(fromCellsCoreInner, knight) 

                upperInnerFromToKnight <- .cs(fromCellsUpperInner, knight[3:8])
                lowerInnerFromToKnight <- .cs(fromCellsLowerInner, knight[1:6])
                leftInnerFromToKnight <- .cs(fromCellsLeftInner, knight[c(1,2,4,6:8)])
                rightInnerFromToKnight <- .cs(fromCellsRightInner, knight[c(1:3,5,7,8)])
                upperleftInnerFromToKnight <- .cs(fromCellUpperleftInner, knight[c(4,6:8)])
                upperrightInnerFromToKnight <- .cs(fromCellUpperrightInner, knight[c(3,5,7,8)])
                lowerleftInnerFromToKnight <- .cs(fromCellLowerleftInner, knight[c(1,2,4,6)])
                lowerrightInnerFromToKnight <- .cs(fromCellLowerrightInner, knight[c(1:3,5)])

                leftOuterFromToKnight <- .cs(fromCellsLeftOuter, knight[c(2,4,6,8)])
                rightOuterFromToKnight <- .cs(fromCellsRightOuter, knight[c(1,3,5,7)])
                upperOuterFromToKnight <- .cs(fromCellsUpperOuter, knight[5:8])
                lowerOuterFromToKnight <- .cs(fromCellsLowerOuter, knight[1:4])
                upperleftUnderFromToKnight <- .cs(fromCellUpperleftUnder, knight[c(4,6,8)])
                upperrightLeftFromToKnight <- .cs(fromCellUpperrightLeft, knight[c(5,7,8)])
                lowerleftUpFromToKnight <- .cs(fromCellLowerleftUp, knight[c(2,4,6)])
                lowerrightUpFromToKnight <- .cs(fromCellLowerrightUp, knight[c(1,3,5)])
                upperleftRightFromToKnight <- .cs(fromCellUpperleftRight, knight[6:8])
                upperrightUnderFromToKnight <- .cs(fromCellUpperrightUnder, knight[c(3,5,7)])
                lowerleftRightFromToKnight <- .cs(fromCellLowerleftRight, knight[c(1,2,4)])
                lowerrightLeftFromToKnight <- .cs(fromCellLowerrightLeft, knight[1:3])
                upperleftFromToKnight <- .cs(fromCellUpperleft, knight[c(6,8)])
                upperrightFromToKnight <- .cs(fromCellUpperright, knight[c(5,7)])
                lowerleftFromToKnight <- .cs(fromCellLowerleft, knight[c(2,4)])
                lowerrightFromToKnight <- .cs(fromCellLowerright, knight[c(1,3)])

                fromto3 <- rbind(coreInnerFromToKnight, upperInnerFromToKnight, lowerInnerFromToKnight, leftInnerFromToKnight, rightInnerFromToKnight, upperleftInnerFromToKnight, upperrightInnerFromToKnight, lowerleftInnerFromToKnight, lowerrightInnerFromToKnight, leftOuterFromToKnight, rightOuterFromToKnight, upperOuterFromToKnight, lowerOuterFromToKnight, upperleftUnderFromToKnight, upperrightLeftFromToKnight, lowerleftUpFromToKnight, lowerrightUpFromToKnight, upperleftRightFromToKnight, upperrightUnderFromToKnight, lowerleftRightFromToKnight, lowerrightLeftFromToKnight, upperleftFromToKnight, upperrightFromToKnight, lowerleftFromToKnight, lowerrightFromToKnight)
                fromto3 <- subset(fromto3,fromto3[,2] %in% toCells)

                if (outerMeridianConnect) {
                        knightLeft <- c(-nCols-1, -2, +2*nCols-2, 3*nCols-1)
                        knightRight <- c(-3*nCols+1, -2*nCols+2, +2, nCols+1)
                        leftInnerFromToKnight <- .cs(fromCellsLeftInner, knightLeft[c(2,3)])
                        rightInnerFromToKnight <- .cs(fromCellsRightInner, knightRight[c(2,3)])
                        upperleftInnerFromToKnight <- .cs(fromCellUpperleftInner, knightLeft[c(2,3)])
                        upperrightInnerFromToKnight <- .cs(fromCellUpperrightInner, knightRight[c(2,3)])
                        lowerleftInnerFromToKnight <- .cs(fromCellLowerleftInner, knightLeft[c(2,3)])
                        lowerrightInnerFromToKnight <- .cs(fromCellLowerrightInner, knightRight[c(2,3)])

                        leftOuterFromToKnight <- .cs(fromCellsLeftOuter, knightLeft)
                        rightOuterFromToKnight <- .cs(fromCellsRightOuter, knightRight)
                        upperleftUnderFromToKnight <- .cs(fromCellUpperleftUnder, knightLeft[2:4])
                        upperrightLeftFromToKnight <- .cs(fromCellUpperrightLeft, knightRight[3])
                        lowerleftUpFromToKnight <- .cs(fromCellLowerleftUp, knightLeft[1:3])
                        lowerrightUpFromToKnight <- .cs(fromCellLowerrightUp, knightRight[1:3])
                        upperleftRightFromToKnight <- .cs(fromCellUpperleftRight, knightLeft[c(3)])
                        upperrightUnderFromToKnight <- .cs(fromCellUpperrightUnder, knightRight[2:4])
                        lowerleftRightFromToKnight <- .cs(fromCellLowerleftRight, knightLeft[2])
                        lowerrightLeftFromToKnight <- .cs(fromCellLowerrightLeft, knightRight[2])
                        upperleftFromToKnight <- .cs(fromCellUpperleft, knightLeft[c(3,4)])
                        upperrightFromToKnight <- .cs(fromCellUpperright, knightRight[c(3,4)])
                        lowerleftFromToKnight <- .cs(fromCellLowerleft, knightLeft[c(1,2)])
                        lowerrightFromToKnight <- .cs(fromCellLowerright, knightRight[c(1,2)])

                        fromto3 <- rbind(fromto3, leftInnerFromToKnight, rightInnerFromToKnight, upperleftInnerFromToKnight, upperrightInnerFromToKnight, lowerleftInnerFromToKnight, lowerrightInnerFromToKnight, leftOuterFromToKnight, rightOuterFromToKnight, upperleftUnderFromToKnight, upperrightLeftFromToKnight, lowerleftUpFromToKnight, lowerrightUpFromToKnight, upperleftRightFromToKnight, upperrightUnderFromToKnight, lowerleftRightFromToKnight, lowerrightLeftFromToKnight, upperleftFromToKnight, upperrightFromToKnight, lowerleftFromToKnight, lowerrightFromToKnight)
                }

                fromto3 <- subset(fromto3,fromto3[,2] %in% toCells)     
                fromto <- rbind(fromto,fromto3)
        }
        colnames(fromto) <- c(from,to)
        return(fromto)
}
.adjBishop <- function(raster, fromCells, toCells, outerMeridianConnect)  {
        nCols <- ncol(raster)
        nCells <- ncell(raster)

        left <- seq(nCols+1,(nCells-2*nCols+1),by=nCols) 
        right <- seq(2*nCols,nCells-nCols,by=nCols)
        upper <- 2:(nCols-1)
        lower <- seq((nCells-nCols+2),(nCells-1),by=1)
        upperleft <- 1
        upperright <- nCols
        lowerleft <- nCells-nCols+1
        lowerright <- nCells
        fromCellsCore <- as.integer(setdiff(fromCells,(c(left,right,upper,lower,upperleft,upperright,lowerleft,lowerright))))
        fromCellsUpper <- as.integer(intersect(fromCells,upper))
        fromCellsLower <- as.integer(intersect(fromCells,lower))
        fromCellsLeft <- as.integer(intersect(fromCells,left))
        fromCellsRight <- as.integer(intersect(fromCells,right))
        fromCellUpperleft <- as.integer(intersect(fromCells,upperleft))
        fromCellUpperright <- as.integer(intersect(fromCells,upperright))
        fromCellLowerleft <- as.integer(intersect(fromCells,lowerleft))
        fromCellLowerright <- as.integer(intersect(fromCells,lowerright))

        bishop <- as.integer(c(-nCols-1, -nCols+1, nCols-1,+nCols+1))

        coreFromToBishop <- .cs(fromCellsCore,bishop)
        upperFromToBishop <- .cs(fromCellsUpper,bishop[3:4])
        lowerFromToBishop <- .cs(fromCellsLower,bishop[1:2])
        leftFromToBishop <- .cs(fromCellsLeft,bishop[c(2,4)])
        rightFromToBishop <- .cs(fromCellsRight,bishop[c(1,3)])
        upperleftFromToBishop <- .cs(fromCellUpperleft,bishop[4])
        upperrightFromToBishop <- .cs(fromCellUpperright,bishop[3])
        lowerleftFromToBishop <- .cs(fromCellLowerleft,bishop[2])
        lowerrightFromToBishop <- .cs(fromCellLowerright,bishop[1])
        fromto <- rbind(coreFromToBishop,upperFromToBishop,lowerFromToBishop,leftFromToBishop,rightFromToBishop,upperleftFromToBishop,upperrightFromToBishop,lowerleftFromToBishop,lowerrightFromToBishop)

        if (outerMeridianConnect) {
                meridianFromLeft <- rbind(
                        .cs(fromCellsLeft,c(2*nCols-1,-1)),
                        cbind(fromCellUpperleft,as.integer(fromCellUpperleft+2*nCols-1)),
                        cbind(fromCellLowerleft,as.integer(fromCellLowerleft-1))
                        ) 
                meridianFromRight <- rbind(
                        cbind(rep(fromCellsRight,times=2),as.integer(c(fromCellsRight-2*nCols+1,fromCellsRight+1))),
                        cbind(fromCellUpperright,as.integer(fromCellUpperright+1)),
                        cbind(fromCellLowerright,as.integer(fromCellLowerright-2*nCols+1))
                        )
                fromto <- rbind(fromto,meridianFromLeft,meridianFromRight)
        }
        fromto <- subset(fromto,fromto[,2] %in% toCells)
        return(fromto)
}

# Author: Robert J. Hijmans
# Date :  September 2011
# Version 1.0
# Licence GPL v3
.adjacentUD <- function(x, cells, ngb, include) {
        # ngb should be a matrix with 
        # one and only one cell with value 0 (the focal cell), 
        # at least one cell with value 1 (the adjacent cells)
        # cells with other values are ignored (not considered adjacent)
        rs <- res(x)
        rn <- raster(ngb)
        center <- which(values(rn)==0)
        if (include) {
                ngb[center] <- 1
        }
        rc <- rowFromCell(rn, center)
        cc <- colFromCell(rn, center)

        xngb <- yngb <- ngb
        xngb[] <- rep(1:ncol(ngb), each=nrow(ngb)) - cc 
        yngb[] <- rep(nrow(ngb):1, ncol(ngb)) - (nrow(ngb)-rc+1)
        ngb[ngb != 1] <- NA
        xngb <- na.omit(as.vector( xngb * rs[1] * ngb))
        yngb <- na.omit(as.vector( yngb * rs[2] * ngb))

        xy <- xyFromCell(x, cells)
        X <- apply(xy[,1,drop=FALSE], 1, function(z) z + xngb )
        Y <- apply(xy[,2,drop=FALSE], 1, function(z) z + yngb )
        c(as.vector(X), as.vector(Y))
}
adjacent <- function(x, cells, directions=4, pairs=TRUE, target=NULL, sorted=FALSE, include=FALSE, id=FALSE) {
        if (is.character(directions)) { 
                directions <- tolower(directions) 
        }
        x <- raster(x)
        r <- res(x)
        xy <- xyFromCell(x, cells)
        mat <- FALSE
        if (is.matrix(directions)) {
                stopifnot(length(which(directions==0)) == 1)
                stopifnot(length(which(directions==1)) > 0)

                d <- .adjacentUD(x, cells, directions, include)

                directions <- sum(directions==1, na.rm=TRUE)
                mat <- TRUE

        } else if (directions==4) {
                if (include) {
                        d <- c(xy[,1], xy[,1]-r[1], xy[,1]+r[1], xy[,1], xy[,1], xy[,2], xy[,2], xy[,2], xy[,2]+r[2], xy[,2]-r[2])
                } else {
                        d <- c(xy[,1]-r[1], xy[,1]+r[1], xy[,1], xy[,1], xy[,2], xy[,2], xy[,2]+r[2], xy[,2]-r[2])
                }

        } else if (directions==8) {
                if (include) {
                        d <- c(xy[,1], rep(xy[,1]-r[1], 3), rep(xy[,1]+r[1],3), xy[,1], xy[,1],
                                 xy[,2], rep(c(xy[,2]+r[2], xy[,2], xy[,2]-r[2]), 2), xy[,2]+r[2], xy[,2]-r[2])
                } else {
                        d <- c(rep(xy[,1]-r[1], 3), rep(xy[,1]+r[1],3), xy[,1], xy[,1],
                                rep(c(xy[,2]+r[2], xy[,2], xy[,2]-r[2]), 2), xy[,2]+r[2], xy[,2]-r[2])
                }
        } else if (directions==16) {
                r2 <- r * 2
                if (include) {
                        d <- c(xy[,1], rep(xy[,1]-r2[1], 2), rep(xy[,1]+r2[1], 2),
                                rep(xy[,1]-r[1], 5), rep(xy[,1]+r[1], 5),
                                xy[,1], xy[,1], 

                                xy[,2], rep(c(xy[,2]+r[2], xy[,2]-r[2]), 2),
                                rep(c(xy[,2]+r2[2], xy[,2]+r[2], xy[,2], xy[,2]-r[2], xy[,2]-r2[2]), 2),
                                xy[,2]+r[2], xy[,2]-r[2])
                } else {
                        d <- c(rep(xy[,1]-r2[1], 2), rep(xy[,1]+r2[1], 2),
                                rep(xy[,1]-r[1], 5), rep(xy[,1]+r[1], 5),
                                xy[,1], xy[,1], 

                                rep(c(xy[,2]+r[2], xy[,2]-r[2]), 2),
                                rep(c(xy[,2]+r2[2], xy[,2]+r[2], xy[,2], xy[,2]-r[2], xy[,2]-r2[2]), 2),
                                xy[,2]+r[2], xy[,2]-r[2])
                }                                       

        } else if (directions=='bishop') {
                if (include) {
                        d <- c(xy[,1], rep(xy[,1]-r[1], 2), rep(xy[,1]+r[1],2), xy[,2], rep(c(xy[,2]+r[2], xy[,2]-r[2]), 2))
                } else {
                        d <- c(rep(xy[,1]-r[1], 2), rep(xy[,1]+r[1],2), rep(c(xy[,2]+r[2], xy[,2]-r[2]), 2))            
                }
                directions <- 4 # to make pairs

        } else {
                stop('directions should be one of: 4, 8, 16, bishop, or a matrix')
        }
        if (include) directions <- directions + 1

        d <- matrix(d, ncol=2)
        if (.isGlobalLonLat(x)) {
                # normalize longitude to -180..180
                d[,1] <- (d[,1] + 180) %% 360 - 180
        }

        if (pairs) {
                if (mat) {
                        cell <- rep(cells, each=directions)             
                } else {
                        cell <- rep(cells, directions)
                }

                if (id) {
                        if (mat) {
                                ID <- rep(1:length(cells), each=directions)
                        } else {
                                ID <- rep(1:length(cells), directions)
                        }
                        d <- na.omit(cbind(ID, cell, cellFromXY(x, d)))
                        attr(d, 'na.action') <- NULL
                        colnames(d) <- c('id', 'from', 'to')
                        if (! is.null(target)) {
                                d <- d[d[,3] %in% target, ]
                        }

                } else {
                        d <- na.omit(cbind(cell, cellFromXY(x, d)))
                        attr(d, 'na.action') <- NULL
                        colnames(d) <- c('from', 'to')
                        if (! is.null(target)) {
                                d <- d[d[,2] %in% target, ]
                        }
                }
                if (sorted) {
                        d <- d[order(d[,1], d[,2]),]
                }
        } else {
                d <- as.vector(unique(na.omit(cellFromXY(x, d))))
                if (! is.null(target)) {
                        d <- intersect(d, target)
                }
                if (sorted) {
                        d <- sort(d)
                }
        }
        d
}

# Author: Robert J. Hijmans
# Date : July 2010
# Version 1.0
# Licence GPL v3
# October 2012: Major overhaul (including C interface)
# November 2012: fixed bug with expand=F
# June 2014: support for aggregation over z (layers) in addition to x and y
setMethod('aggregate', signature(x='Raster'), 
function(x, fact=2, fun='mean', expand=TRUE, na.rm=TRUE, filename=, ...)  {
        doC <- list(...)$doC
        if (is.null(doC)) {
                doC <- TRUE
        }
        nl <- nlayers(x)
        fact <- round(fact)
        lf <- length(fact)
        if (lf == 1) {
                fact <- c(fact, fact, 1)
        } else if (lf == 2) {
                fact <- c(fact, 1)
        } else if (lf > 3) {
                stop('fact should have length 1, 2, or 3')
        }
        if (nl < 2) {
                fact[3] <- 1
        }
        if (any(fact < 1)) {
                stop('fact should be > 0')
        }
        if (! any(fact > 1)) {
                stop('fact should be > 1')
        }
        xfact <- fact[1]
        yfact <- fact[2]
        zfact <- fact[3]

        ncx <- ncol(x)
        nrx <- nrow(x)
        if (xfact > ncx) {
                warning('aggregation factor is larger than the number of columns') 
                xfact <- ncx
        }
        if (yfact > nrx) {
                warning('aggregation factor is larger than the number of rows')
                yfact <- nrx
        }
        if (zfact > nl) {
                warning('aggregation factor is larger than the number of layers')
                zfact <- nl
        }
        addlyrs <- 0
        if (expand) {
                rsteps <- as.integer(ceiling(nrx/yfact))
                csteps <- as.integer(ceiling(ncx/xfact))
                lsteps <- as.integer(ceiling(nl/zfact))

                lastcol <- ncx
                lastrow <- nrx
                lastlyr <- lsteps * zfact
                if (lastlyr > nl ) {
                        addlyrs <- lastlyr - nl
                }
                lyrs <- 1:nl

                #addcols <- csteps * xfact - ncx
                #addrows <- rsteps * yfact - nrx
        } else  {
                rsteps <- as.integer(floor(nrx/yfact))
                csteps <- as.integer(floor(ncx/xfact))
                lsteps <- as.integer(floor(nl/zfact))

                lastcol <- min(csteps * xfact, ncx)
                lastrow <- min(rsteps * yfact, nrx)
                lastlyr <- min(lsteps * zfact, nl)

                lyrs <- 1:lastlyr
        }

        ymn <- ymax(x) - rsteps * yfact * yres(x)
        xmx <- xmin(x) + csteps * xfact * xres(x)

        if (lsteps > 1) {
                out <- brick(x, values=FALSE)
        } else {
                out <- raster(x)                
        }
        extent(out) <- extent(xmin(x), xmx, ymn, ymax(x))
        dim(out) <- c(rsteps, csteps, lsteps) 
        ncout <- ncol(out)
        if (zfact == 1) {
                names(out) <- names(x)
        }


        if (! hasValues(x) ) {  
                return(out) 
        }       
        fun <- .makeTextFun(fun)
        if (class(fun) == 'character') { 
                op <- as.integer(match(fun, c('sum', 'mean', 'min', 'max')) - 1)
        } else {
                op <- NA
        }

        if (zfact > 1) {
                xyzfact <- xfact*yfact*zfact
                dims <- as.integer(c(lastrow, lastcol, nl+addlyrs, xfact, yfact, zfact))
                if ( canProcessInMemory(x)) {
                        v <- getValuesBlock(x, 1, lastrow, 1, lastcol, lyrs)
                        if (addlyrs > 0) {
                                add <- matrix(NA, nrow=nrow(v), ncol=addlyrs)
                                v <- cbind(v, add)
                        }
                        v <- .Call(aggregate_get, as.double(v), as.integer(dims), PACKAGE='raster')
                        v <- matrix(v, nrow=xyzfact)
                        v <- apply(v, 2, fun, na.rm=na.rm)
                        out <- setValues(out, v)
                        if (filename != '') {
                                out <- writeRaster(out, filename, ...)
                        }
                        return(out)
                } else {
                        xx <- brick(x, values=FALSE)
                        if (!expand) {
                                nrow(xx) <- (nrow(x) %/% yfact) * yfact
                        }               
                        tr <- blockSize(xx, minrows=yfact)
                        st <- round(tr$nrows[1] / yfact) * yfact
                        tr$n <- ceiling(lastrow / st)
                        tr$row <- c(1, cumsum(rep(st, tr$n-1))+1)
                        tr$nrows <- rep(st, tr$n)
                        tr$write <- cumsum(c(1, ceiling(tr$nrows[1:(tr$n-1)]/yfact)))
                        tr$nrows[tr$n] <-  nrow(xx) - tr$row[tr$n] + 1
                        tr$outrows <- ceiling(tr$nrows/yfact)

                        pb <- pbCreate(tr$n, label='aggregate', ...)
                        x <- readStart(x, ...)  
                        out <- writeStart(out, filename=filename, ...)
                        for (i in 1:tr$n) {
                                dims[1] <- as.integer(tr$nrows[i])
                                vals <- getValuesBlock(x, tr$row[i], tr$nrows[i], 1, lastcol, lyrs)
                                if (addlyrs > 0) {
                                        add <- rep(NA, nrow(vals)*addlyrs)
                                        vals <- c(vals, add)
                                }
                                vals <- .Call(aggregate_get, as.double(vals), as.integer(dims), PACKAGE='raster')
                                vals <- matrix(vals, nrow=xyzfact)
                                vals <- apply(vals, 2, fun, na.rm=na.rm)
                                out <- writeValues(out, matrix(vals, ncol=nl), tr$write[i])
                                pbStep(pb, i) 
                        }
                        pbClose(pb)
                        out <- writeStop(out)
                        x <- readStop(x)
                        return(out)     
                }
        }

        if (!is.na(op) & doC) {

                if ( canProcessInMemory(x)) {

                        dims <- as.integer(c(lastrow, lastcol, nl, dim(out)[1:2], xfact, yfact))
                        x <- getValuesBlock(x, 1, lastrow, 1, lastcol)
                        out <- setValues(out, .Call(aggregate, as.double(x), op, as.integer(na.rm), dims, PACKAGE='raster'))
                        if (filename != '') {
                                out <- writeRaster(out, filename, ...)
                        }
                        return(out)

                } else {

                        xx <- brick(x, values=FALSE)
                        if (!expand) {
                                nrow(xx) <- (nrow(x) %/% yfact) * yfact
                        }               
                        tr <- blockSize(xx, minrows=yfact)
                        st <- round(tr$nrows[1] / yfact) * yfact
                        tr$n <- ceiling(lastrow / st)
                        tr$row <- c(1, cumsum(rep(st, tr$n-1))+1)
                        tr$nrows <- rep(st, tr$n)
                        tr$write <- cumsum(c(1, ceiling(tr$nrows[1:(tr$n-1)]/yfact)))
                        tr$nrows[tr$n] <-  nrow(xx) - tr$row[tr$n] + 1
                        tr$outrows <- ceiling(tr$nrows/yfact)

                        pb <- pbCreate(tr$n, label='aggregate', ...)
                        x <- readStart(x, ...)  
                        dims <- as.integer(c(lastrow, lastcol, nl, dim(out)[1:2], xfact, yfact))
                        out <- writeStart(out, filename=filename, ...)
                        if (inherits(out, 'RasterBrick')) {
                                for (i in 1:tr$n) {
                                        dims[c(1, 4)] = as.integer(c(tr$nrows[i], tr$outrows[i]))
                                        vals <- getValuesBlock(x, tr$row[i], tr$nrows[i], 1, lastcol)
                                        vals <- .Call(aggregate, as.double(vals), op, as.integer(na.rm), dims, PACKAGE='raster')
                                        out <- writeValues(out, matrix(vals, ncol=nl), tr$write[i])
                                        pbStep(pb, i) 
                                }
                        } else {
                                for (i in 1:tr$n) {
                                        dims[c(1, 4)] = as.integer(c(tr$nrows[i], tr$outrows[i]))
                                        vals <- getValuesBlock(x, tr$row[i], tr$nrows[i], 1, lastcol)
                                        vals <- .Call(aggregate, as.double(vals), op,   as.integer(na.rm), dims, PACKAGE='raster')
                                        out <- writeValues(out, vals, tr$write[i])
                                        pbStep(pb, i) 
                                }
                        }
                        pbClose(pb)
                        out <- writeStop(out)
                        x <- readStop(x)
                        return(out)     
                }
        }

 # else not implemented in C  

        if (nl < 2) {   
                if (class(fun) == 'character') { 
                        rowcalc <- TRUE 
                        fun <- .getColFun(fun)
                } else { 
                        rowcalc <- FALSE 
                }

                if ( canProcessInMemory(x)) {
                        if (expand) {
                                m <- ceiling(nrx / yfact)
                        } else {
                                m <- floor(nrx / yfact)
                        }
                        vv <- matrix(NA, nrow= yfact*xfact, ncol=csteps * m)
                        vend <- 0
                        vvstart <- 1

                        if (expand) {
                                vals <- getValues(x)
                                yf <- nrx %% yfact
                        } else {
                                vals <- getValuesBlock(x, 1, lastrow, 1, lastcol)
                                yf <- 0
                        }
                        for (j in 1:m) {
                                if (j == m & yf > 0) {
                                        vstart <- vend + 1
                                        vend <- vend + (lastcol * yf)
                                        mv <- matrix(vals[vstart:vend], nrow=yf, byrow=TRUE )
                                        temp <- matrix(nrow=yf*xfact, ncol=csteps)
                                        temp[1:length(mv)] <- mv
                                        cols <- 1:(csteps) + (m-1) * csteps
                                        vv[1:nrow(temp), cols] <- temp

                                } else {
                                        vstart <- vend + 1
                                        vend <- vend + (lastcol * yfact)
                                        mv <- matrix(vals[vstart:vend], nrow=yfact, byrow=TRUE )

                                        vv[vvstart:(vvstart+length(mv)-1)] <- as.vector(mv)
                                        vvstart <- vvstart + ncout*nrow(vv)
                                }
                        }
                        if (rowcalc) {
                                vals <- fun(vv, na.rm=na.rm )
                        } else {
                                vals <- apply(vv, 2, fun, na.rm=na.rm )
                        }
                        out <- setValues(out, as.vector(vals))
                        if (filename != '') {
                                out <- writeRaster(out, filename, ...)
                        }
                        return(out)

                } else {

                        out <- writeStart(out, filename=filename, ...)


                        tr <- blockSize(x, minrows=yfact)
                        st <- round(tr$nrows[1] / yfact) * yfact
                        tr$n <- ceiling(lastrow / st)
                        tr$row <- c(1, cumsum(rep(st, tr$n-1))+1)
                        tr$nrows <- rep(st, tr$n)
                        tr$write <- cumsum(c(1, ceiling(tr$nrows[1:(tr$n-1)]/yfact)))
                        dif <- sum(tr$nrows) - nrow(x)
                        if (dif > 0) {
                                if (expand) {
                                        tr$nrows[tr$n] <-  tr$nrows[tr$n] - dif
                                } else {
                                        dif <- dif %/% xfact
                                        if (dif > 0) {
                                                tr$nrows[tr$n] <- dif * xfact
                                        } else {
                                                tr$n <- tr$n - 1
                                        }
                                }
                        }

                        pb <- pbCreate(tr$n, label='aggregate', ...)
                        x <- readStart(x, ...)  
                        m <- tr$nrows[1] / yfact
                        vv <- matrix(NA, nrow= yfact*xfact, ncol=csteps * m)

                        w <- getOption('warn')
                        on.exit(options('warn' = w))
                        options('warn'=-1) 

                        for (i in 1:(tr$n-1)) {
                                vals <- getValuesBlock(x, tr$row[i], tr$nrows[i], 1, lastcol)           
                                vend <- 0
                                vvstart <- 1
                                for (j in 1:m) {
                                        vstart <- vend + 1
                                        vend <- vend + (lastcol * yfact)
                                        mv <- matrix(vals[vstart:vend], nrow=yfact, byrow=TRUE )

                                        vv[vvstart:(vvstart+length(mv)-1)] <- as.vector(mv)                     
                                        vvstart <- vvstart + ncout*nrow(vv)
                                }
                                if (rowcalc) {
                                        vals <- fun(vv, na.rm=na.rm )
                                } else {
                                        vals <- apply(vv, 2, fun, na.rm=na.rm )
                                }
                                out <- writeValues(out, vals, tr$write[i])
                                pbStep(pb, i) 
                        } 
        #       if (i==tr$n) { 
                        i <- tr$n
                        vals <- getValuesBlock(x, tr$row[i], tr$nrows[i], 1, lastcol)           
                        m <- ceiling(tr$nrows[i] / yfact)
                        vv <- matrix(NA, nrow= yfact*xfact, ncol=csteps * m)
                        vend <- 0
                        vvstart <- 1
                        yf <- tr$nrows[i] %% yfact
                        for (j in 1:m) {
                                if (j == m & yf > 0) {
                                        vstart <- vend + 1
                                        vend <- vend + (lastcol * yf)
                                        mv <- matrix(vals[vstart:vend], nrow=yf, byrow=TRUE )
                                        temp <- matrix(nrow=yf*xfact, ncol=csteps)
                                        temp[1:length(mv)] <- mv
                                        cols <- 1:(csteps) + (m-1) * csteps
                                        vv[1:nrow(temp), cols] <- temp

                                } else {
                                        vstart <- vend + 1
                                        vend <- vend + (lastcol * yfact)
                                        mv <- matrix(vals[vstart:vend], nrow=yfact, byrow=TRUE )

                                        vv[vvstart:(vvstart+length(mv)-1)] <- as.vector(mv)
                                        vvstart <- vvstart + ncout*nrow(vv)
                                }
                        }
                        if (rowcalc) {
                                vals <- fun(vv, na.rm=na.rm )
                        } else {
                                vals <- apply(vv, 2, fun, na.rm=na.rm )
                        }
                        pbStep(pb, i) 
                        out <- writeValues(out, vals, tr$write[i])
                        pbClose(pb)
                        x <- readStop(x)        
                        out <- writeStop(out)
                        return(out)
                }

        } else { # nlayers > 1

                if (canProcessInMemory(x, nlayers(x)+2)) {

                        if (class(fun) == 'character') { 
                                op <- as.integer(match(fun, c('sum', 'mean', 'min', 'max')) - 1)
                        }
                        if (!is.na(op) & doC) {
                                dim <- c(dim(x), dim(out)[1:2], xfact, yfact)
                                v  <- .Call(aggregate, 
                                                as.double(getValues(x)), op, as.integer(na.rm), 
                                                as.integer(dim), PACKAGE='raster')

                                out <- setValues(out, matrix(v, ncol=dim[3]))           
                                return(out)     
                        }

                        xx <- raster(x)         
                        x <- getValues(x)
                        cols <- rep(rep(1:csteps, each=xfact)[1:ncol(xx)], times=nrow(xx))
                        rows <- rep(1:rsteps, each=ncol(xx) * yfact)[1:ncell(xx)]
                        cells <- cellFromRowCol(xx, rows, cols)

                        x <- as.matrix( aggregate(x, list(cells), fun, na.rm=na.rm ))[,-1]
                        rm(cells)

                        x <- setValues(out, x)
                        if (filename != ) {
                                x <- writeRaster(x, filename=filename, ...)
                        }
                        return(x)
                } else  { 

                        cols <- rep(rep(1:csteps,each=xfact)[1:ncol(x)], times=yfact)
                        rows <- rep(1, each=(ncol(x) * yfact))

                        out <- writeStart(out, filename=filename, ...)
                        x <- readStart(x, ...)  

                        cells <- cellFromRowCol(x, rows, cols)
                        nrows <- yfact
                        w <- getOption('warn')
                        on.exit(options('warn' = w))
                        options('warn'=-1) 

                        pb <- pbCreate(rsteps, label='aggregate', ...)
                        for (r in 1:rsteps) {
                                startrow <- 1 + (r - 1) * yfact
                                if ( r==rsteps) {
                                        endrow <- min(nrow(x), startrow + yfact - 1)
                                        nrows <- endrow - startrow + 1
                                        theserows <- (startrow * rows)[1:(ncol(x)*nrows)]
                                        cols <- cols[1:(ncol(x)*nrows)]
                                        cells <- cellFromRowCol(x, theserows, cols)
                                }       
                                vals <- getValues(x, startrow, nrows)
                                vals <- as.matrix( aggregate(vals, list(cells), fun, na.rm=na.rm ))[,-1]

                                out <- writeValues(out, vals, r)
                                pbStep(pb, r) 
                        } 
                        pbClose(pb)
                        out <- writeStop(out)
                        x <- readStop(x)        
                        return(out)
                }       
        }
}
)
#library(raster)
#r <- raster(nc=9, nr=9)
#r <- raster()
#r[] = 1:ncell(r)
#.aggtest(r, 5, 'min', doC=T)
#aggregate(s, c(2,1,3), 'min', expand=F)

# Author: Robert J. Hijmans
# Date : July 2010
# Version 1.0
# Licence GPL v3
.aggregate_old <- function(x, fact=2, fun=mean, expand=TRUE, na.rm=TRUE, filename=, ...)  {
        if (length(fact)==1) {
                fact <- as.integer(round(fact))
                if (fact < 2) { stop('fact should be > 1') }
                xfact <- yfact <- fact
        } else if (length(fact)==2) {
                xfact <- as.integer(round(fact[[1]]))
                yfact <- as.integer(round(fact[[2]]))
                if (xfact < 2) { stop('fact[[1]] should be > 1') } 
                if (yfact < 2) { stop('fact[[2]] should be > 1') }
        } else {
                stop('length(fact) should be 1 or 2')
        }
        if (xfact > ncol(x)) {
                warning('aggregation factor is larger than the number of columns') 
                xfact <- ncol(x)
        }
        if (yfact > nrow(x)) {
                warning('aggregation factor is larger than the number of rows')
                yfact <- nrow(x)
        }
        if (expand) {
                rsteps <- as.integer(ceiling(nrow(x)/yfact))
                csteps <- as.integer(ceiling(ncol(x)/xfact))
                lastcol <- x@ncols
                lastrow <- x@nrows
        } else  {
                rsteps <- as.integer(floor(nrow(x)/yfact))
                csteps <- as.integer(floor(ncol(x)/xfact))
                lastcol <- min(csteps * xfact, x@ncols)
                lastrow <- min(rsteps * yfact, x@nrows)
        }

        ymn <- ymax(x) - rsteps * yfact * yres(x)
        xmx <- xmin(x) + csteps * xfact * xres(x)

        nl <- nlayers(x)
        if (nl > 1) {
                out <- brick(x, values=FALSE)
        } else {
                out <- raster(x)                
        }
        extent(out) <- extent(xmin(x), xmx, ymn, ymax(x))
        dim(out) <- c(rsteps, csteps) 
        names(out) <- names(x)
        if (! hasValues(x) ) {  return(out) }   

        if (nl < 2) {   
                fun <- .makeTextFun(fun)
                if (class(fun) == 'character') { 
                        rowcalc <- TRUE 
                        fun <- .getRowFun(fun)
                } else { 
                        rowcalc <- FALSE 
                }
                if (! canProcessInMemory(x)) {
                        if (filename == '') { 
                                filename <- rasterTmpFile() 
                        }
                }

                if (filename == '') {
                        v <- matrix(NA, ncol=nrow(out), nrow=ncol(out))
                } else {
                        out <- writeStart(out, filename=filename, ...)
                }

                pb <- pbCreate(rsteps, ...)

                        #vv <- matrix(ncol= csteps * yfact, nrow=rsteps * xfact)
                vv <- matrix(nrow= yfact * xfact, ncol=csteps)
                w <- getOption('warn')
                on.exit(options('warn' = w))
                options('warn'=-1) 
                for (r in 1:rsteps) {

                        startrow <- 1 + (r - 1) * yfact
                        vals <- getValuesBlock(x, startrow, yfact, 1, lastcol)

                        if (r==rsteps) { 
                                endrow <- min(x@nrows, (startrow + yfact - 1))
                                nrows <- endrow - startrow + 1
                                vals <- matrix(vals, nrow=nrows, byrow=TRUE )
                                vv[] <- NA 
                                vvv <- vv[1:(nrows*xfact), ,drop=FALSE]
                                vvv[1:length(vals)] <- vals
                                vv[1:nrow(vvv),] <- vvv

                        } else {
                                vals <- matrix(vals, nrow=yfact, byrow=TRUE )
                                vv[1:length(vals)] = vals
                        }

                        if (rowcalc) {
                                vals <- fun(t(vv), na.rm=na.rm )
                        } else {
                                vals <- apply(vv, 2, fun, na.rm=na.rm )
                        }
                        if (filename == ) {
                                v[, r] <- vals
                        } else {
                                out <- writeValues(out, vals, r)
                        }
                        pbStep(pb, r) 
                } 

                pbClose(pb)
                if (filename == ) { 
                        values(out) <- as.vector(v)
                } else {
                        out <- writeStop(out)
                }
                return(out)

        } else { # nlayers > 1


                if (canProcessInMemory(x, nlayers(x)+2)) {

                        xx <- raster(x)         
                        x <- getValues(x)
                        cols <- rep(rep(1:csteps, each=xfact)[1:ncol(xx)], times=nrow(xx))
                        rows <- rep(1:rsteps, each=ncol(xx) * yfact)[1:ncell(xx)]
                        cells <- cellFromRowCol(xx, rows, cols)

                        x <- as.matrix( aggregate(x, list(cells), fun, na.rm=na.rm ))[,-1]
                        rm(cells)

                        x <- setValues(out, x)
                        if (filename != ) {
                                x <- writeRaster(x, filename=filename, ...)
                        }
                        return(x)
                } else  { 

                        cols <- rep(rep(1:csteps,each=xfact)[1:ncol(x)], times=yfact)
                        rows <- rep(1, each=(ncol(x) * yfact))
                        cells <- cellFromRowCol(x, rows, cols)
                        nrows <- yfact
                        w <- getOption('warn')
                        on.exit(options('warn' = w))
                        options('warn'=-1) 

                        out <- writeStart(out, filename=filename, ...)
                        pb <- pbCreate(rsteps, ...)
                        for (r in 1:rsteps) {
                                startrow <- 1 + (r - 1) * yfact
                                if ( r==rsteps) {
                                        endrow <- min(nrow(x), startrow + yfact - 1)
                                        nrows <- endrow - startrow + 1
                                        theserows <- (startrow * rows)[1:(ncol(x)*nrows)]
                                        cols <- cols[1:(ncol(x)*nrows)]
                                        cells <- cellFromRowCol(x, theserows, cols)
                                }       
                                vals <- getValues(x, startrow, nrows)
                                vals <- as.matrix( aggregate(vals, list(cells), fun, na.rm=na.rm ))[,-1]

                                out <- writeValues(out, vals, r)
                                pbStep(pb, r) 
                        } 
                        pbClose(pb)
                        out <- writeStop(out)
                        return(out)
                }

        }
}

# Author: Robert J. Hijmans
# Date : July 2010
# Version 1.0
# Licence GPL v3
# October 2012: Major overhaul (including C interface)
# November 2012: fixed bug with expand=F
setMethod('aggregate', signature(x='Raster'), 
function(x, fact=2, fun='mean', expand=TRUE, na.rm=TRUE, filename=, ...)  {
        doC <- list(...)$doC
        if (is.null(doC)) {
                doC <- TRUE
        }
        fact <- rep(as.integer(round(fact)), length.out=2)
        xfact <- fact[1]
        yfact <- fact[2]
        if (xfact < 1 | yfact < 1) { stop('fact should be > 0') } 
        if (xfact < 2 & yfact < 2) { stop('fact[1] or fact[2] should be > 1') } 

        if (xfact > ncol(x)) {
                warning('aggregation factor is larger than the number of columns') 
                xfact <- ncol(x)
        }
        if (yfact > nrow(x)) {
                warning('aggregation factor is larger than the number of rows')
                yfact <- nrow(x)
        }
        ncx <- ncol(x)
        nrx <- nrow(x)
        if (expand) {
                rsteps <- as.integer(ceiling(nrx/yfact))
                csteps <- as.integer(ceiling(ncx/xfact))
                lastcol <- x@ncols
                lastrow <- x@nrows
                #addcols <- csteps * xfact - ncx
                #addrows <- rsteps * yfact - nrx
        } else  {
                rsteps <- as.integer(floor(nrx/yfact))
                csteps <- as.integer(floor(ncx/xfact))
                lastcol <- min(csteps * xfact, x@ncols)
                lastrow <- min(rsteps * yfact, x@nrows)
        }


        ymn <- ymax(x) - rsteps * yfact * yres(x)
        xmx <- xmin(x) + csteps * xfact * xres(x)

        nl <- nlayers(x)
        if (nl > 1) {
                out <- brick(x, values=FALSE)
        } else {
                out <- raster(x)                
        }
        extent(out) <- extent(xmin(x), xmx, ymn, ymax(x))
        dim(out) <- c(rsteps, csteps) 
        names(out) <- names(x)
        ncout <- ncol(out)
        if (! hasValues(x) ) {  return(out) }   
        fun <- .makeTextFun(fun)
        if (class(fun) == 'character') { 
                op <- as.integer(match(fun, c('sum', 'mean', 'min', 'max')) - 1)
        } else {
                op <- NA
        }

        if (!is.na(op) & doC) {

                if ( canProcessInMemory(x)) {

                        dims <- as.integer(c(lastrow, lastcol, nl, dim(out)[1:2], xfact, yfact))
                        x <- getValuesBlock(x, 1, lastrow, 1, lastcol)
                        out <- setValues(out, .Call(aggregate, as.double(x), op, as.integer(na.rm), dims, PACKAGE='raster'))
                        if (filename != '') {
                                out <- writeRaster(out, filename, ...)
                        }
                        return(out)

                } else {

                        xx <- brick(x, values=FALSE)
                        if (!expand) {
                                xx <- brick(x, values=FALSE)
                                nrow(xx) <- (nrow(x) %/% yfact) * yfact
                        }               
                        tr <- blockSize(xx, minrows=yfact)
                        st <- round(tr$nrows[1] / yfact) * yfact
                        tr$n <- ceiling(lastrow / st)
                        tr$row <- c(1, cumsum(rep(st, tr$n-1))+1)
                        tr$nrows <- rep(st, tr$n)
                        tr$write <- cumsum(c(1, ceiling(tr$nrows[1:(tr$n-1)]/yfact)))
                        tr$nrows[tr$n] <-  nrow(xx) - tr$row[tr$n] + 1
                        tr$outrows <- ceiling(tr$nrows/yfact)

                        pb <- pbCreate(tr$n, label='aggregate', ...)
                        x <- readStart(x, ...)  
                        dims <- as.integer(c(lastrow, lastcol, nl, dim(out)[1:2], xfact, yfact))
                        out <- writeStart(out, filename=filename, ...)
                        if (inherits(out, 'RasterBrick')) {
                                for (i in 1:tr$n) {
                                        dims[c(1, 4)] = as.integer(c(tr$nrows[i], tr$outrows[i]))
                                        vals <- getValuesBlock(x, tr$row[i], tr$nrows[i], 1, lastcol)
                                        vals <- .Call(aggregate, as.double(vals), op, as.integer(na.rm), dims, PACKAGE='raster')
                                        out <- writeValues(out, matrix(vals, ncol=nl), tr$write[i])
                                        pbStep(pb, i) 
                                }
                        } else {
                                for (i in 1:tr$n) {
                                        dims[c(1, 4)] = as.integer(c(tr$nrows[i], tr$outrows[i]))
                                        vals <- getValuesBlock(x, tr$row[i], tr$nrows[i], 1, lastcol)
                                        vals <- .Call(aggregate, as.double(vals), op,   as.integer(na.rm), dims, PACKAGE='raster')
                                        out <- writeValues(out, vals, tr$write[i])
                                        pbStep(pb, i) 
                                }
                        }
                        pbClose(pb)
                        out <- writeStop(out)
                        x <- readStop(x)
                        return(out)     
                }
        }

 # else not implemented in C  

        if (nl < 2) {   
                if (class(fun) == 'character') { 
                        rowcalc <- TRUE 
                        fun <- .getColFun(fun)
                } else { 
                        rowcalc <- FALSE 
                }

                if ( canProcessInMemory(x)) {
                        if (expand) {
                                m <- ceiling(nrx / yfact)
                        } else {
                                m <- floor(nrx / yfact)
                        }
                        vv <- matrix(NA, nrow= yfact*xfact, ncol=csteps * m)
                        vend <- 0
                        vvstart <- 1

                        if (expand) {
                                vals <- getValues(x)
                                yf <- nrx %% yfact
                        } else {
                                vals <- getValuesBlock(x, 1, lastrow, 1, lastcol)
                                yf <- 0
                        }
                        for (j in 1:m) {
                                if (j == m & yf > 0) {
                                        vstart <- vend + 1
                                        vend <- vend + (lastcol * yf)
                                        mv <- matrix(vals[vstart:vend], nrow=yf, byrow=TRUE )
                                        temp <- matrix(nrow=yf*xfact, ncol=csteps)
                                        temp[1:length(mv)] <- mv
                                        cols <- 1:(csteps) + (m-1) * csteps
                                        vv[1:nrow(temp), cols] <- temp

                                } else {
                                        vstart <- vend + 1
                                        vend <- vend + (lastcol * yfact)
                                        mv <- matrix(vals[vstart:vend], nrow=yfact, byrow=TRUE )

                                        vv[vvstart:(vvstart+length(mv)-1)] <- as.vector(mv)
                                        vvstart <- vvstart + ncout*nrow(vv)
                                }
                        }
                        if (rowcalc) {
                                vals <- fun(vv, na.rm=na.rm )
                        } else {
                                vals <- apply(vv, 2, fun, na.rm=na.rm )
                        }
                        out <- setValues(out, as.vector(vals))
                        if (filename != '') {
                                out <- writeRaster(out, filename, ...)
                        }
                        return(out)

                } else {

                        out <- writeStart(out, filename=filename, ...)


                        tr <- blockSize(x, minrows=yfact)
                        st <- round(tr$nrows[1] / yfact) * yfact
                        tr$n <- ceiling(lastrow / st)
                        tr$row <- c(1, cumsum(rep(st, tr$n-1))+1)
                        tr$nrows <- rep(st, tr$n)
                        tr$write <- cumsum(c(1, ceiling(tr$nrows[1:(tr$n-1)]/yfact)))
                        dif <- sum(tr$nrows) - nrow(x)
                        if (dif > 0) {
                                if (expand) {
                                        tr$nrows[tr$n] <-  tr$nrows[tr$n] - dif
                                } else {
                                        dif <- dif %/% xfact
                                        if (dif > 0) {
                                                tr$nrows[tr$n] <- dif * xfact
                                        } else {
                                                tr$n <- tr$n - 1
                                        }
                                }
                        }

                        pb <- pbCreate(tr$n, label='aggregate', ...)
                        x <- readStart(x, ...)  
                        m <- tr$nrows[1] / yfact
                        vv <- matrix(NA, nrow= yfact*xfact, ncol=csteps * m)

                        w <- getOption('warn')
                        on.exit(options('warn' = w))
                        options('warn'=-1) 

                        for (i in 1:(tr$n-1)) {
                                vals <- getValuesBlock(x, tr$row[i], tr$nrows[i], 1, lastcol)           
                                vend <- 0
                                vvstart <- 1
                                for (j in 1:m) {
                                        vstart <- vend + 1
                                        vend <- vend + (lastcol * yfact)
                                        mv <- matrix(vals[vstart:vend], nrow=yfact, byrow=TRUE )

                                        vv[vvstart:(vvstart+length(mv)-1)] <- as.vector(mv)                     
                                        vvstart <- vvstart + ncout*nrow(vv)
                                }
                                if (rowcalc) {
                                        vals <- fun(vv, na.rm=na.rm )
                                } else {
                                        vals <- apply(vv, 2, fun, na.rm=na.rm )
                                }
                                out <- writeValues(out, vals, tr$write[i])
                                pbStep(pb, i) 
                        } 
        #       if (i==tr$n) { 
                        i <- tr$n
                        vals <- getValuesBlock(x, tr$row[i], tr$nrows[i], 1, lastcol)           
                        m <- ceiling(tr$nrows[i] / yfact)
                        vv <- matrix(NA, nrow= yfact*xfact, ncol=csteps * m)
                        vend <- 0
                        vvstart <- 1
                        yf <- tr$nrows[i] %% yfact
                        for (j in 1:m) {
                                if (j == m & yf > 0) {
                                        vstart <- vend + 1
                                        vend <- vend + (lastcol * yf)
                                        mv <- matrix(vals[vstart:vend], nrow=yf, byrow=TRUE )
                                        temp <- matrix(nrow=yf*xfact, ncol=csteps)
                                        temp[1:length(mv)] <- mv
                                        cols <- 1:(csteps) + (m-1) * csteps
                                        vv[1:nrow(temp), cols] <- temp

                                } else {
                                        vstart <- vend + 1
                                        vend <- vend + (lastcol * yfact)
                                        mv <- matrix(vals[vstart:vend], nrow=yfact, byrow=TRUE )

                                        vv[vvstart:(vvstart+length(mv)-1)] <- as.vector(mv)
                                        vvstart <- vvstart + ncout*nrow(vv)
                                }
                        }
                        if (rowcalc) {
                                vals <- fun(vv, na.rm=na.rm )
                        } else {
                                vals <- apply(vv, 2, fun, na.rm=na.rm )
                        }
                        pbStep(pb, i) 
                        out <- writeValues(out, vals, tr$write[i])
                        pbClose(pb)
                        x <- readStop(x)        
                        out <- writeStop(out)
                        return(out)
                }

        } else { # nlayers > 1

                if (canProcessInMemory(x, nlayers(x)+2)) {

                        if (class(fun) == 'character') { 
                                op <- as.integer(match(fun, c('sum', 'mean', 'min', 'max')) - 1)
                        }
                        if (!is.na(op) & doC) {
                                dim <- c(dim(x), dim(out)[1:2], xfact, yfact)
                                v  <- .Call(aggregate, 
                                                as.double(getValues(x)), op, as.integer(na.rm), 
                                                as.integer(dim), PACKAGE='raster')

                                out <- setValues(out, matrix(v, ncol=dim[3]))           
                                return(out)     
                        }

                        xx <- raster(x)         
                        x <- getValues(x)
                        cols <- rep(rep(1:csteps, each=xfact)[1:ncol(xx)], times=nrow(xx))
                        rows <- rep(1:rsteps, each=ncol(xx) * yfact)[1:ncell(xx)]
                        cells <- cellFromRowCol(xx, rows, cols)

                        x <- as.matrix( aggregate(x, list(cells), fun, na.rm=na.rm ))[,-1]
                        rm(cells)

                        x <- setValues(out, x)
                        if (filename != ) {
                                x <- writeRaster(x, filename=filename, ...)
                        }
                        return(x)
                } else  { 

                        cols <- rep(rep(1:csteps,each=xfact)[1:ncol(x)], times=yfact)
                        rows <- rep(1, each=(ncol(x) * yfact))

                        out <- writeStart(out, filename=filename, ...)
                        x <- readStart(x, ...)  

                        cells <- cellFromRowCol(x, rows, cols)
                        nrows <- yfact
                        w <- getOption('warn')
                        on.exit(options('warn' = w))
                        options('warn'=-1) 

                        pb <- pbCreate(rsteps, label='aggregate', ...)
                        for (r in 1:rsteps) {
                                startrow <- 1 + (r - 1) * yfact
                                if ( r==rsteps) {
                                        endrow <- min(nrow(x), startrow + yfact - 1)
                                        nrows <- endrow - startrow + 1
                                        theserows <- (startrow * rows)[1:(ncol(x)*nrows)]
                                        cols <- cols[1:(ncol(x)*nrows)]
                                        cells <- cellFromRowCol(x, theserows, cols)
                                }       
                                vals <- getValues(x, startrow, nrows)
                                vals <- as.matrix( aggregate(vals, list(cells), fun, na.rm=na.rm ))[,-1]

                                out <- writeValues(out, vals, r)
                                pbStep(pb, r) 
                        } 
                        pbClose(pb)
                        out <- writeStop(out)
                        x <- readStop(x)        
                        return(out)
                }       
        }
}
)
#library(raster)
#r <- raster(nc=9, nr=9)
#r <- raster()
#r[] = 1:ncell(r)
#.aggtest(r, 5, 'min', doC=T)